利用多个PPAC的x,y信息进行traking的原理如图所示,图中红点代表探测器测量的位置。由于PPAC的探测效率,对于每个入射粒子只有部分PPAC的x或y方向可给出位置信息。假设对于某一个入射粒子,有两个以上(含两个)探测器给出x位置信息 (如图:1Ax,1Bx,3x) 时,入射粒子的x-z径迹可由上述x-z点进行线性拟合得到。y-z径迹的处理方法与x-z径迹的处理方法一致。利用上述线性方程,可内插或外推某一z位置的x和y信息。图中空心黑点代表由线性方程计算得到的位置。
![]()
1) 若F8PPAC3中有有效的位置而且F8PPAC1和F8PPAC2中至少有一层PPAC有有效的位置。
2) 若F8PPAC3没有有效的位置,则


通过准直孔将粒子准直,入射到探测器上。此时测得的位置分辨率$\sigma$与探测器本征位置分辨率$\sigma_0$的关系为 $$ \sigma^2=\sigma^2_0+\sigma^2_{geo} $$ ,其中$\sigma_{geo}$为准直孔引入的分辨率(图a)。
上述方法一般用放射源进行,而PPAC探测器的位置分辨与入射粒子的种类以及能量相关,因此由放射源得到的分辨率不能直接用在束流条件下。



假设有$N$个粒子穿过探测器灵敏面积,探测器实际探测到的数目为$N_{det}$, 则探测效率$\epsilon$为 $$ \epsilon = \frac{N_{det}}{N} $$
在本实验设置下,可以有很多方法求探测器的效率:
例1.ppac的阳极信号的计数作为分母,阴极计数作为分子(图a)。
例2.选择任意两组(i,j)探测器,拟合得到径迹, 记录径迹穿过探测器k灵敏面积的所有事件数目$N_{ij}$。记录在上述条件下探测器k具有正确位置信息的事件数目$N_{kij}$(图b)。 $$ \epsilon = \frac{N_{kij}}{N_{ij}} $$
用上述做法可分别计算探测器的x,y,x-y的效率
下面示例代码假设同时有1A,2A,3探测器给出x,y的位置信息。
$PPACF8[i][j]$ -Calibrated Data(位置已刻度好)
| Branch | PPAC | code |
|---|---|---|
| PPACF8[0][0] | PPAC 1 Layer A X (mm) | xx[0] |
| PPACF8[0][1] | PPAC 1 Layer A Y (mm) | yy[0] |
| PPACF8[0][2] | PPAC 1 Layer A Z for X-plan (mm) | xz[0] |
| PPACF8[0][3] | PPAC 1 Layer A Z for Y-plan (mm) | yz[0] |
| PPACF8[0][4] | PPAC 1 Layer A Anode time (ns) | |
| PPACF8[1][0] | PPAC 1 Layer B X (mm) | |
| PPACF8[1][1] | PPAC 1 Layer B Y (mm) | |
| PPACF8[1][2] | PPAC 1 Layer B Z for X-plan (mm) | |
| PPACF8[1][3] | PPAC 1 Layer B Z for Y-plan (mm) | |
| PPACF8[1][4] | PPAC 1 Layer B Anode time (ns) | |
| PPACF8[2][0-4] | PPAC 2 Layer A * | xx[1],yy[1],xz[1],yz[1] |
| PPACF8[3][0-4] | PPAC 2 Layer B * | xx2b,yy2b,xz2b,yz2b |
| PPACF8[4][0-4] | PPAC 3 * | xx[2],yy[2],xz[2],yz[2] |
root -l f8ppac001.root
>tree->MakeClass("tracking");
>.q
编辑 tracking.h和tracking.C 保存
使用方法:
root -l
> .L tracking.C
> tracking t
>t.Loop()
>.q
root -l tracking.root //分析
在traking.h 代码中增加了用户变量,成员函数SetBranch(),TrackInit() 以及SetTrace()定义
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Wed Mar 11 09:58:39 2020 by ROOT version 6.18/04
// from TTree tree/tree
// found on file: f8ppac001.root
//////////////////////////////////////////////////////////
#ifndef tracking_h
#define tracking_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
// Header file for the classes stored in the TTree if any.
class tracking {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
Int_t fCurrent; //!current Tree number in a TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Float_t PPACF8[5][5];
Float_t F8PPACRawData[5][5];
Int_t beamTrig;
Int_t must2Trig;
Float_t targetX,targetY;
//by user
Double_t xx[3],xz[3],yy[3],yz[3];//1A,2A,3
Double_t xx2b[2],yy2b[2],xz2b,yz2b;//2B x,y, 0-measured, 1- fitted.
Double_t dx[3],dy[3];//residual
Double_t tx,ty;//target position
Double_t c2nx,c2ny;//chi2/ndf for xfit,yfit
// List of branches
TBranch *b_PPACF8; //!
TBranch *b_F8PPACRawData; //!
TBranch *b_beamTrig; //!
TBranch *b_must2Trig; //!
TBranch *b_targetX; //!
TBranch *b_targetY; //!
tracking(TTree *tree=0);
virtual ~tracking();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
virtual Long64_t LoadTree(Long64_t entry);
virtual void Init(TTree *tree);
virtual void Loop();
virtual void SetBranch(TTree *tree);//by user
virtual void TrackInit();//by user
virtual void SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max);//by user
virtual Bool_t Notify();
virtual void Show(Long64_t entry = -1);
};
#endif
#ifdef tracking_cxx
tracking::tracking(TTree *tree) : fChain(0)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
if (tree == 0) {
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("f8ppac001.root");
if (!f || !f->IsOpen()) {
f = new TFile("f8ppac001.root");
}
f->GetObject("tree",tree);
}
Init(tree);
}
tracking::~tracking()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t tracking::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t tracking::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->GetTreeNumber() != fCurrent) {
fCurrent = fChain->GetTreeNumber();
Notify();
}
return centry;
}
void tracking::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("PPACF8", PPACF8, &b_PPACF8);
fChain->SetBranchAddress("F8PPACRawData", F8PPACRawData, &b_F8PPACRawData);
fChain->SetBranchAddress("beamTrig", &beamTrig, &b_beamTrig);
fChain->SetBranchAddress("must2Trig", &must2Trig, &b_must2Trig);
fChain->SetBranchAddress("targetX",&targetX,&b_targetX);
fChain->SetBranchAddress("targetY",&targetY,&b_targetY);
Notify();
}
Bool_t tracking::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
void tracking::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
Int_t tracking::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
#endif // #ifdef tracking_cxx
在tracking的基础上进行了修改
#define tracking_cxx
#include "tracking.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TFitResult.h>
void tracking::SetBranch(TTree *tree)
{
//measured pos
tree->Branch("xx",&xx,"xx[3]/D");//1A,2A,3
tree->Branch("xz",&xz,"xz[3]/D");
tree->Branch("yy",&yy,"yy[3]/D");
tree->Branch("yz",&yz,"yz[3]/D");
//difference between measured and calculated -for pos resolution.
tree->Branch("dx",&dx,"dx[3]/D");
tree->Branch("dy",&dy,"dy[3]/D");
//2B x,y
tree->Branch("xx2b",&xx2b,"xx2b[2]/D");
tree->Branch("yy2b",&yy2b,"yy2b[2]/D");
//target x-y
tree->Branch("tx",&tx,"tx/D");
tree->Branch("ty",&ty,"ty/D");
//ch2/ndf for linear fitting.
tree->Branch("c2nx",&c2nx,"c2nx/D");
tree->Branch("c2ny",&c2ny,"c2ny/D");
tree->Branch("beamTrig",&beamTrig,"beamTrig/I");
tree->Branch("must2Trig",&must2Trig,"must2Trig/I");
tree->Branch("targetX",&targetX,"targetX");
tree->Branch("targetY",&targetY,"targetY");
}
void tracking::TrackInit()
{
tx=-999;
ty=-999;
//1A
xx[0]=PPACF8[0][0];
yy[0]=PPACF8[0][1];
xz[0]=PPACF8[0][2];
yz[0]=PPACF8[0][3];
//2A
xx[1]=PPACF8[2][0];
yy[1]=PPACF8[2][1];
xz[1]=PPACF8[2][2];
yz[1]=PPACF8[2][3];
//3
xx[2]=PPACF8[4][0];
yy[2]=PPACF8[4][1];
xz[2]=PPACF8[4][2];
yz[2]=PPACF8[4][3];
//2B
xx2b[0]=PPACF8[3][0];
yy2b[0]=PPACF8[3][1];
xz2b=PPACF8[3][2];
yz2b=PPACF8[3][3];
xx2b[1]=-1000;
yy2b[1]=-1000;
}
void tracking::SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max){
if(h==0) return;
if(min>=max) return;
for(int i=min;i<max;i++){
h->Fill(i,(Int_t)(i*k+b));
}
}
void tracking::Loop()
{
TH2D *htf8xz=new TH2D("htf8xz","xz trace by ppac",2200,-2000,200,300,-150,150);
TH2D* htf8yz=new TH2D("htf8yz","yz trace by ppac",2200,-2000,200,300,-150,150);
TFile *opf=new TFile("tracking.root","recreate");
TTree *tree=new TTree("tree","ppac traking");
SetBranch(tree);
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
TrackInit();
bool b1a=abs(xx[0])<150 && abs(yy[0])<150;
bool b2a=abs(xx[1])<150 && abs(yy[1])<150;
bool b3=abs(xx[2])<150 && abs(yy[2])<150;
if(!b1a || !b2a || !b3) continue;
//fit x-z trajectory
TFitResultPtr r;
TGraph *grx=new TGraph(3,xz,xx);
TF1 *fx=new TF1("fx","pol1",-2000,0);
r=grx->Fit(fx,"SQ");
xx2b[1]=fx->Eval(xz2b);
tx=fx->Eval(0);
SetTrace(htf8xz,fx->GetParameter(1),fx->GetParameter(0),-1800,0);
for(int i=0;i<3;i++) dx[i]=xx[i]-fx->Eval(xz[i]);
c2nx=r->Chi2()/r->Ndf();
delete grx;
delete fx;
//fit y-z trajectory
TGraph *gry=new TGraph(3,yz,yy);
TF1 *fy=new TF1("fy","pol1",-2000,0);
r=gry->Fit(fy,"SQ");
yy2b[1]=fy->Eval(yz2b);
ty=fy->Eval(0);
SetTrace(htf8yz,fy->GetParameter(1),fy->GetParameter(0),-1800,0);
for(int i=0;i<3;i++) dy[i]=yy[i]-fy->Eval(yz[i]);
c2ny=r->Chi2()/r->Ndf(); //对于任何程序的自动拟合,原则上都要输出拟合误差进行评估
delete gry;
delete fy;
tree->Fill();
if(jentry%10000==0) cout<<"processing "<<jentry<<endl;
}
htf8xz->Write();
htf8yz->Write();
tree->Write();
opf->Close();
}
%jsroot on
TFile *ipf=new TFile("tracking.root");
TTree *tree=(TTree*) ipf->Get("tree");
TCanvas *c1=new TCanvas("c1","c1");
TH2D *hxz=(TH2D*) ipf->Get("htf8xz");
hxz->Draw("colz");
c1->Draw();
TH2D *hyz=(TH2D*) ipf->Get("htf8yz");
hyz->Draw("colz");
c1->Draw();
tree->Draw("ty:tx>>htx(120,-60,60,120,-60,60)","must2Trig","colz");
c1->Draw();
tree->Draw("tx:ty>>htx(120,-60,60,120,-60,60)","beamTrig","colz");
c1->Draw();
$\chi^2/Ndf : tx$
tree->Draw("dx[0]>>(200,-5,5)");
c1->Draw();//residual
tree->Draw("dx[0]:dx[1]");
c1->Draw();
tree->Draw("c2ny:dy[0]>>hh(40,-10,10,200,0,1000)","","colz");
c1->Draw();//从chi2/ndf图上可看出,部分事件的径迹拟合误差很大,这部分要在后续数据处理中去掉。
tree->Draw("c2ny>>hh(200,0,1000)","","");
c1->Draw();//从chi2/ndf图上可看出,部分事件的径迹拟合误差很大,这部分要在后续数据处理中去掉。
tree->Draw("tx:ty>>(120,-60,60)","c2nx<10 && c2ny<10 && beamTrig ","colz");
c1->Draw();//
tree->Draw("tx:ty>>(120,-60,60,120,-60,60)","(c2nx>20 || c2ny>20) && beamTrig ","colz");
c1->Draw();//
Long64_t N_track;
Long64_t Nx_det, Ny_det,Nxy_det;
TCut c2btrack="abs(xx2b[1])<100 && abs(yy2b[1])<100";//径迹穿过探测器的灵敏面积
TCut c2bx="abs(xx2b[0])<100";//x面有正确信号
TCut c2by="abs(yy2b[0])<100";//y面有正确信号
tree->Draw("yy2b[1]:xx2b[1]>>(200,-100,100,200,-100,100)",c2btrack,"colz");//fitted
N_track=tree->GetEntries(c2btrack);//得到给定条件下的计数。
cout<<N_track<<endl;
c1->Draw();
tree->Draw("yy2b[0]:xx2b[0]>>(200,-100,100,200,-100,100)",c2bx&&c2by&&c2btrack,"colz");
Nx_det=tree->GetEntries(c2bx && c2btrack);
Ny_det=tree->GetEntries(c2by && c2btrack);
Nxy_det=tree->GetEntries(c2bx && c2by && c2btrack);
cout<<Nx_det<<" "<<Ny_det<<" "<<Nxy_det<<endl;
c1->Draw();
Double_t ex,ey,exy;
ex=Double_t(Nx_det)/N_track;
ey=Double_t(Ny_det)/N_track;
exy=Double_t(Nxy_det)/N_track;
TString eff;
eff.Form("PPAC2B:\n eff_x=%.2f%%, \n eff_y=%.2f%%,\n eff_x*eff_y=%.2f%%, \n eff_xy=%.2f%%",ex*100,ey*100,ex*ey*100,exy*100);
cout<<eff.Data()<<endl;